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ABSTRACT 

We attempt to verify recent claims (made using semi-analytic models) that for the collapse of 
spherical homogeneous molecular clouds, fragmentation of the self-gravitating disc that sub- 
sequently forms can be predicted using the cloud's initial angular momentum alone. In effect, 
this condition is equivalent to requiring the resulting disc be sufficiently extended in order 
to fragment, in line with studies of isolated discs. We use smoothed particle hydrodynamics 
with hybrid radiative transfer to investigate this claim, confirming that in general, homoge- 
neous spherical molecular clouds will produce fragmenting self-gravitating discs if the ratio 
of rotational kinetic energy to gravitational potential energy is greater than fa 5 x 10^"^, where 
this result is relatively insensitive to the initial thermal energy. This condition begins to fail 
at higher cloud masses, suggesting that sufficient mass at large radii governs fragmentation. 
While these results are based on highly idealised initial conditions, and may not hold if the 
disc's accretion from the surrounding envelope is sufficiently asymmetric, or if the density 
structure is perturbed, they provide a sensible lower limit for the minimum angular momen- 
tum required to fragment a disc in the absence of significant external turbulence. 

Key words: stars: formation, accretion, accretion discs, methods: numerical, radiative trans- 
fer, hydrodynamics 



1 INTRODUCTION 

\ Cold, dense molecular cloud cores are though t to be the sites 
. of low-mass star formation l iTerebev et al. l [T984h . Typically, these 
■ cores will contain an excess of angular momentum (c ompared with 
the rotational angular momentum of low-mass stars, cf lCaselli et al.l 
[2002). If these cores are sufficiently dense, they will collapse under 
their own gravity to form a protostar plus protostellar disc, with a 
substantial envelope surrounding both. In this pre-main sequence 
phase, we can expect that most of the star's eventual mass will be 
processed by the protostellar disc (which in turn accretes this mass 
from a surrounding envelope). 

The star's accretion is inextricably linked with the pro- 
cess of outward angular momentum transport - understanding the 
mechanisms by which angular momentum transport functions in 
these early phases is therefore essential to theories of star for- 
mation. While turbulence induced by the magneto-rotational in- 
stability (MRI) define s angular momentum transp ort at late times 
( iBlaes & Balbu3ll994lB"albus & Papaloizoulll999h . the disc is un- 
likely to be even weakly ionised at early times. We must therefore 
seek another means by which to generate angular momentum trans- 
port. 

Simulations suggest that the disc is likely to be self- 
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gravitating at these early ti mes, a s the disc mass will be large 
relati ve to the star mass teatd |2010| ; iMachida & Matsumotd 
l201lh . If the disc becomes gravitationally unstable, this may pro- 
vide the dominant transport me chanism th r ough the growth of 
self-regulated gravito-turbulence jBosjl 19841 : iLin & Pringldl 19871 : 
iLaughlin & Bodenheimeil [l994h . A self-gravitating disc will be- 
come subject to gravitational instability (GI) if the Toomre param- 
eter, Q ( Toomre 1964): 

(1) 



^ ttGE 

where Cs is the sound speed, E is the surface density, and k is the 
epicyclic frequency (in the case of a Keplerian disc, k is equal to 
the angular frequency fi). If Q < 1, axisymmetric perturbations 
will grow within the disc. Numerical simulations have established 
that discs wit h Q < 1.5— 1. 7 will be unstable to non-axisymmetric 
perturbations jPurisen et al . 2007). The onset of GI is accompanied 
by the growth of quasi-steady spira l structures, which d rive shocks 
with Mach numbers of order unity jCossins et al ]|2009l) . 

The growth of GI in a self-gravitating disc will gener- 
ally lead to a state commonly referred to as marginal instabil- 
ity jPaczvnskH [79781 where the heating due to GI is balanced 
by radiative cooling. The combination of these two processes 
produces a quasi-equilibrium state where the instability's am- 
plitudes are modulated, providing a self-regulated outward an- 
gular momentum transport process ( IGammi 3 l200ll : |pickettetal.l 
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I2OO3I ; iLodato & Ricell2004l ; ICai et alj|2008h . The growth and de- 
cay of these structures provides the so-called gravito-turbulence, 
allowing the evolutio n of the disc to be described in a pseudo- 
viscous fashion (e.g. Balbus & Papaloizou 1999; Lodato & Ricd 
|2004 iRice & Armitagell2009l : Iciarkell2009h . where the transport 
can be character ised using the dimension less turbulent viscos- 
ity parameter a dShakura & Sunvaev|[l973h . While this local a- 
parametrisation breaks down in certain limits iForgan et al.ll201 ih . 
it remains a useful tool to compare to other transport mechanisms 
such as the MRI. 

Self-gravitating disc evolution is also an important compo- 
nent of planet formation theory, not merely th rough its effects 
on the growth and accretion of planetary cores jRice et ai] |2004 
IClarke & Lodata2009,) . If the disc is self-gravitating, its dynamical 
evolution can r apidly form g i ant planets and brown dwarfs by dis c 
fragmentation ( lKuipeij[l95ll : IBossIHwtI : IStamatellos et al.ll2007h . 
While the Q criterion is a necessary condition for disc fragmenta- 
tion, it is insufficient. A second criterion is required, relating to the 
local cooling time: 

(2) 



du/dt' 



where u is the specific internal energy. It is often more suitable to 
discuss cooling in terms of a dimensionless parameter, fie '■ 

Pc = fcoolfi. (3) 

iGammid lliioih showed using shearing-sheet simulations that if 
/3c = const., fragmentation can only occur if jSc < /3crit, where 
/3crit is approximately 3, although it was later shown that /3crit de- 
pends on the local equati on of state, in pa rticular the value of the 
ratio of specific heats, 7 iRice et aT]|2005l) . The correct values of 
/?crit(7 ) have been challenged by recent work, and ma y be slightly 
higher jMeni & Bat"3l201 1 j ; |PaardekooDer etal]|201 ll) . In the case 
of smoothed particle hydrodynamics (SPH) simulations as are car- 
ried out in this work, this appears to be due to resolution effects 
pushing the convergen ce limit to higher particle number than was 
previously anticipated jLodato & Clarke 2011), as well as possible 
issues regarding the correct implementation of cooling prescrip- 
tions (Rice, Forgan and Armitage, submitted). 

Alternatively, the second fragmentation criterion can be 
couched in terms of the local stresses in the disc. In the Shakura- 
Sunyaev formalism, the effective turbulent viscosity v = aCsH, 
where H is the disc scale height. Stresses have two principal origins 
in self-gravitating discsQ- Reynolds stresses (induced by correlated 
velocity perturbations) and gravitational stresses (induced by cor- 
related perturbations in the gravitational potential). Assuming that 
the disc is in thermal equilibrium, the value of a can be determined 
by balancing the pseudo- viscous heating and the radiative cooling: 



dlnfl 



(4) 



dlnr J 7(7 — l)/3c ' 

iRice et al show that a single critical value of q = Ocrit ~ 

0.06 can be used as the second fragmentation criterion, for any 
value of 7. If the disc is irradiated, it has been shown that 
a, does not uniquely determines the likelihood of fragmenta- 
tion, and that short cooling ti mes can induce fragmentation even 
in low-stress environments ( (Rice et al.ll201 1[) . Studies have how- 
ever shown that this Ocrit value may not apply if /3c can vary 



dClarke et al.l |2007|) , or if the disc's accretion from the enve- 
lope is substantial dKratter et alj I2OO8I : iHarsono et all l201lh . In 
these circumstances, more sophisticated criteria are required (see 
e.g. Iforgan & Rice 201 1.) . In spite of these circumstances, it has 
been established by many authors that the requisite condi tions for 
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2009; Clarke 2009; Rice et al. 2010). These results hold in semi- 
analytic calculations, grid-based simulations and particle-based 
simulations. This does not limit the usefulness of disc frag- 
mentation as an explanation for observable exoplanets - indeed, 
it explains the properties of systems such as HR 8799 when 
core accretion theories s truggle to do so jKratteret alj I2OI0I ; 
iNero & BiorkmanI l2009t) . iMeru & Bate! ( 1201 id) presented SPH 
simulations which displayed inner disc fragmentation for steep sur- 
face density profiles, but it has been shown that these can be ex- 
plained by failu re to address resolution criteria relating to the arti- 
ficial viscosity iLodato &C1^|20 lib - indeed, a failure that all 
previous SPH simulations of this type have shared. 

As has been stated, the self-gravitating phase is expected to oc- 
cur at early times, when the disc is still embedded in the envelope. 
Therefore, to fully understand the fragmentation of self-gravitating 
discs, it is necessary to include the enve lope's contribution, both in 
its ability to supply the disc with matter dKratter et al.l2008l) . and in 
the resulting stresses it can induce at the disc's surface and beyond 
(Harsono et al. 2011). These stresses result in extensive non-local 
angular momentum transport, a phenomenon whic h will tend to in- 
validate semi-analytic models dForgan et al.ll201 fh . 

A straightforward means by which the envelope can be in- 
cluded is to begin with a rotating molecular cloud as initial condi- 
tions. Frequently, the cloud parameters (j) and A are usec0 where: 
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where iJthorm, E^ot and E^i^v are the cloud's total thermal, ro- 
tational and gravitational energies respectively. For an initially 
isothermal cloud in rigid body rotation, these parameters essentially 
define the cloud's initial temperature and angular velocity: 
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where iici is the cloud's initial radius, k is Boltzmann's constant, ^ 
is the mean mol ecular weight o f the gas and mn is the mass of the 
hydrogen atom. iTohlind dl98ll) considers the collapse of adiabatic 
spheroids (P = Kp^), and defines a series of criteria for these 
clouds to become unstable to nonaxisymmetric perturbations, de- 
pendent on the polytropic index F of the cloud: 
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(9) 



^ A third source of disc stress is magnetic fields, but as the self-gravitating 
phase is typically veiy weakly ionised, the magnitude of these stresses is 
small, and hence we do not consider them here 



^ In the literature, these parameters are known as a and /3 respectively. We 
relabel them to avoid confusion with the subsequent disc parameters. 



Cloud Angular Momentum and Disc Fragmentation 3 



In reality, a collapsing molecular cloud will have a more com- 
plex thermal history that requires a more detailed understanding 
of the eq uation of st ate as a function of the local density and tem- 
perature^ ^ce^^n( 201 G), using a model initially developed by 



lLin&Pringlel l ll99(]|) . carried out semi-analytic calculations of uni- 
form, spherical molecular cloud collapse, with a more complex 
equation of state, using the pseudo-viscous approximation to evolve 
the resulting disc, which continues to accrete from its envelope. 
They show that the resulting self-gravitating discs are subject to 
fragmentation (i.e. they contain regions where acooi > 0.06) when 
A> 1x10"^, and never fragment when A ^ 1 x 10 regardless 
of the initial thermal energy of the cloud. 

It is not clear whether these results can be expected to hold, 
given that non-local angular momentum transport (which must play 
a role, especially at early times) is not accounted for, as well as po- 
tentially asymmetric accretion flows. In addition, these initial con- 
ditions are extremely simple, and not likely to be reproduced in 
reality - a more appropri ate choice w ould be supercritical Bonnor- 
Ebert sp heres (e.g.,Walch et alj2009l) . as well as the addition of tur- 
bulence dWalch et al.ll20IOl) . Finally, as has been previously stated, 
using acrit = 0.06 is probably not appropriate for accreting discs. 

We aim to perform a set of numerical experiments using 
smoothed particle hydrodynamics with hybrid radiative transfer 
jForgan et alj2009l) to appr aise the val i dity o f the ang ular momen- 
tum c ondition set forth by iRice et al] i2010l) . Unlike IWalchetal.1 
(|2009^. we do not hold 4> constant, and therefore we can also de- 
tect whether any dependence on thermal energy exists. Also, we 
can model non-local angular momentum transport, and can mea- 
sure fragmentation directly rather than relying on uncertain criteria, 
allowing us to test more reliably whether an angular momentum 
criterion truly exists. 



2 METHOD 

2.1 SPH and the Hybrid Radiative Transfer Approximation 

Smoothed Particle Hydr o dynamics (SPH ) fcucvl 1 1 9771 : 
iGingold & MonaghanI 1 1971 iMonaghanI 1 19921) is a Lagrangian 
formalism that represents a fluid by a distribution of particles. Each 
particle is assigned a mass, position, internal energy and velocity. 
State variables such as density a nd pressure are then ca lculated 
by interpolation (see reviews by iMonaghanI [l99i, |2005|) . In the 
simulations presented here, the gas is modelled using 500,000 
SPH particles while the star is represented by a point mass particle 
onto which gas particles can accrete, if the y are sufficiently close 
(in our case, within 1 au) and are bound ( lBateetalj[T99a) . An 
SPH particle is converted to a pointmass if the local density has 
exceeded 10~^^ gcm~^, and 

(i) There exists one Jeans mass within the particle's neighbour 
sphere 

(ii) The neighbour sphere is bound (i.e. the total energy inside 
the sphere is negative) 

(iii) The local gravitational potential energy exceeds both the 
thermal and rotational kinetic energy 

The SPH cod e used in this w ork is based on the SPH 
code developed by Bate et af] jl995 !b which uses individual par- 
ticle timesteps, and individually variable smoothing lengths, hi, 
such that the number of nearest neighbours for each particle is 
50 ± 20. The code uses a hybrid method of approximate ra- 
diative transfer jForganetal]|2009l) . which is built on two pre- 



existing radiative algorithms: the polytropic cooling approxima- 
tion d evised by[s tamatellos et al. (20 oj), and flux- lim ited diffusion 
(e.g., Whitehous e & Batei2004l : lMaver et al.l2007 '. see'porgan et alj 
2009 for details). This union allows the effects of both global cool- 
ing and radiative transport to be modelled, without imposing extra 
boundary conditions. 

The opacity and temperature of the gas is calculated using a 
non-trivial equation of state. This accounts for the effects of H2 dis- 
sociation, H" ionisation. He" and He+ ionisation, ice evaporation, 
dust sublimation, molecular absorp tion, bound-free and free-free 
transitions and electron scattering teell&Lirj [19941 :1 Bolev et al.l 
2007; Stam atellos eTai]|2007h . Heating of the disc is also achieved 
by P dV work and shock heating. 



2.2 Initial Cloud Conditions 

The molecular clouds were initialised with 500, 000 SPH particles. 
While this ma y seem under-resolved in the light of the results from 
[Meru & Bati j2011ah . it is important to note that they use a funda- 
mentally different cooling method where particle cooling rates are 
not affected by their neighbours. The hybrid method, through its 
flux limited diffusion component, does account for the neighbour 
spher e. A more appropriate resolution criterion (Lodato & Clarl^ 
I2OIII based on the results of iMeru & Batel I2OI lah shows that 
500,000 particles is sufficient to resolve fragmentation, and is dis- 
cussed in the following section. The particles were distributed in 
a sphere with uniform density. The clouds were initialised with 
varying values of cloud mass Md and A - the cloud radii Rd 
were specified in order to maintain the same initial density po = 
1.46 X 10~^^ gcm~^. This ensures that all clouds have the same 
free-fall time 
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(10) 



The initial cloud temperature is fixed at 5 K, equal to the back- 
ground temperature (beyond which the cloud cannot cool). Note by 
fixing the parameters in this fashion, we allow (j) to vary. We choose 
this setup a s it is similar to that of the semi-analytic calculations of 
iRice et al.l l l201oh . who also allow to vary in this fashion. Their 
results appear to be independent of 0, and we wish to assess if this 
is true in our case. We do not apply a density pertur bation to th e 
clouds, as is often done in simulations of this type (e.g. lBossll986l) . 
The amplitude of the m = 2 mode initially is zero in our simula- 
tions - the subsequent growth of m = 2 spiral modes is purely a 
result of disc torques forming from gravitational instabilities. 



2.3 Resolution Conditions 

To ensure that potential fragmentation is resolved, the minimum 
resolvable Jeans mass - one neighbour group of SPH particles 
l lBate&Burker3l 19971) . around 50 in the case of the code used - 
must be sufficiently small: 



Mmin = 2iVnoi„hmi = 2Mt< 
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The minimum resolvable Jeans mass ranges between SSMg for 
Mci = 0.5A'/q to 120M® for M^i = 2A'/q. As it is expected 
that fragment masses wil l be typically several orders of magnitud e 
higher than these values jKratter et al.l201(il : lForgan &Ricell201lb . 
this establishes that the simulations would comfortably resolve disc 
fragmentation if it were to occur. 
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Table 1. Summary of the cloud parameters investigated in this work. 



Simulation 


Mci (M©) 


Rcl (au) 


A 


</> 


1 


1.0 


2133.33 


4 X lO"-^ 


0.2 


2 


1.0 


2133.33 


5 xlO-3 


0.2 


3 


1.0 


2133.33 


6 xlO-3 


0.2 


4 


1.0 


2133.33 


7 xlO-3 


0.2 


5 


1.0 


2133.33 


8 xlO-3 


0.2 


6 


1.0 


2133.33 


9 xl0"3 


0.2 


7 


0.5 


1693.23 


4 xlO-3 


0.4 


g 


0.5 


1693.23 


5 xlO-3 


0.4 


9 


0.5 


1693.23 


6 xlO-3 


0.4 


10 


0.5 


1693.23 


7 xlO-3 


0.4 


11 


0.5 


1693.23 


8 xlQ--^ 


0.4 


12 


0.5 


1693.23 


9 xlO"^ 


0.4 


13 


2.0 


2687.83 


4 xlO-3 


0.1 


14 


2.0 


2687.83 


5 xlO-3 


0.1 


15 


2.0 


2687.83 


6 xlO-3 


0.1 


16 


2.0 


2687.83 


7 xlO-3 


0.1 


17 


2.0 


2687.83 


8 xlO-3 


0.1 


18 


2.0 


2687.83 


9 xlQ-^ 


0.1 



Perhaps mor e important are the reso lution issues raised by ar- 
tificial viscosity ( lLodato&Clarkell20T ij). While required by the 
SPH code, this artificial viscosity must be quantified so that we 
know where in the disc the artificial viscosity is likely to be lower 
than the effective viscosity generated by the gravitational instabili- 
ties. If the stress induced by artificial viscosity exceeds that of the 
stresses induced by the gravitational instability in som e region of 
the di sc, then fragmentation is likely to be suppressed l lwhitworthi 
Il998l) . 

The linear term for the artificial viscos i ty can be expressed as 
jArtvmowicz & Lubowlll994 lMurravllT99^ : iLodato & Pricel201(]h 



fart = —OtSPHCsh, 



(12) 



where Cs is the local sound speed, h is the local SPH smoothing 
length, and qsph is the linear viscosity coefficient used by the 
SPH code (taken to be 0.1). We can define an effective a param- 
eter associated wi th the artificial viscosity jLodato & R"ic3l2004 
iForgan et al]l201lh 



(13) 



and hence combining eq u ations iT2l a n d l|13t gives 
dArtymowicz & Lubowl 1 19941 : iMurravl 1 19961 : iLodato & Pric3 

mm 



(14) 



This shows that where the vertical structure is not well resolved 
(i.e., is large), artificial viscosity will dominate. In the simula- 
tions presented here, this is likely to be the case inside ~ 10 au, 
so any data inside this region can not be used. We therefore only 
consider results outside this radius, and therefore cannot comment 
on works which demonstrate fragmentation inside this radius. In 
the regions where fragmentation is found, the artificial viscosity is 
less than 1% of th e total, and therefore sati sfies the "less than 5%" 
criterion given by ILodato & Clarkd ( l20Ilh for correctly resolving 
fragmentation. 

Finally, we must recognise that these simulations contain Pois- 
son noise in the initial conditions. After approximately one free-fall 
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Figure 1. The fragmentation boundary in the simulations conducted. The 
boundary between disc fragmentation and no fragmentation depends only 
weakly on initial cloud mass. 



time, the star-disc system begins to settle into a marginally unsta- 
ble state, effectively erasing the initial conditions. Fragmentation 
due to disc instability typically takes an extra 0.5f// after disc for- 
mation, corresponding to many outer rotation periods (ORPs), so 
we do not expect this initial noise to significantly alter our result. 
To confirm this, all simulations were repeated, using a different set 
of randomly generated initial conditions, to confirm that this initial 
noise does not affect the final result. 



3 RESULTS 

3.1 General Trends 

Each of the clouds goes through a set of generic evo- 
lutionary phases, similar i n nature to those described by 
iMasunaga & Inutsukal j2000l) for a non-rotating homogeneous 
cloud. Initially, the collapse is isothermal - despite the central den- 
sity increasing by several orders of magnitude, the central tempera- 
ture remains constant at the background temperature of 5K. As the 
central density increases past 10~^^ gcm~^, the gas becomes op- 
tically thick, and the central temperature begins to rise. As the tem- 
perature rises, rotational degrees of freedom in H2 are activated, 
and the contraction slows somewhat as the cloud develops signif- 
icant thermal pressure support, forming the "first core" (cf lLarsoii 
li969l). As these simulations contain a sufficient amount of initial 
angular momentum, this oblate first core will subs equently form 
the protostellar disc jMachida & Matsumotolbol ih . This process 
begins with the formati on of a bar within the core, which decays 
into spiral structures (cf lBatell2O10h . The second core forms at the 
centre of this disc structure at approximately 2000 K, when H2 
begins to dissociate, removing thermal energy and hence pressure 
support. At this stage, the simulation forms a pointmass. The disc- 
to-pointmass mass ratio is initially larger than 3, and the star ac- 
cretes rapidly. After this rapid accretion phase, the mass ratio drops 
below unity, and the marginally unstable self-gravitating phase be- 
gins (except in the simulations that produce multiple discs). 
Figure [T] displays the end state of all the simulations run in this 
work. The end states can be categorised into one of three broad 
classes: 

(i) Non-fragmenting, marginally unstable self-gravitating discs, 
(denoted by in Figure[TJ, 
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(ii) Self-gravitating discs which fragment several ORPs after 
formation (denoted by triangles in Figure^, 

(iii) Systems which fragment into multiple systems before a sin- 
gle disc system can form (denoted by squares in Figure[TJ. 

It can be seen that the delineation of the three classes shows only 
a very weak dependence on the initial cloud mass Md- Note that 
all clouds have the same initial free-fall time, and the same external 
background temperature. The clouds are initially isothermal, and 
will therefore cool quickly to the background temperature, remov- 
ing any dependence o n initial therma l energy (unlike the adiabatic 
clouds of for example iTohlind Tl 98 ll) The fact that initial angular 
momentum provides a (relatively) simple criterion 

A > 5 X 10"^ (15) 

(and its subsequent weak violation at larger cloud masses) suggests 
that fragmentation requires placing a sufficient quantity of material 
at large radii at early times, having the resulting effect of reducing 
the Q parameter without triggering a strong increase in the local 
cooling time. This allows sufficient stress to be generated as a re- 
sult of the gravitational instability, which fragments the disc. Our 
results contrast with those from simulatio ns where an azimuthal 
density perturbation is applied, for example lBossI ( Il986l) . who find 
that binary formation occurs at A = 4 x 10~^. This discrepancy 
seems reasonable again under the interpretation of sufficient mass 
at large radii at early times promoting fragmentation. 

Discs which do not fragment are qualitatively different in their 
properties to discs that do. These differences are generic to all ini- 
tial cloud masses. For the sake of brevity, we will only explore the 
Mci = IMq clouds in detail. 



3.2 The A/ci = IM© Simulations 

3.2.1 General Properties 

Figure |2] shows surface density plots for Simulations 1, 2, 3, 4, 5 
and 6, after at least 2t// of evolution, or during fragmentation in 
the case of Simulations 3, 4, 5 and 6. In Simulations 3, 4 and 5, 
fragmentation occurs around t — 1.5i//, in Simulation 6 it occurs 
around t = It ff .It is clear that increasing the initial cloud angular 
momentum increases the overall disc extent. There is also a general 
tendency for fragments to form closer to the central object as the 
angular momentum is increased. In all cases, it is apparent that the 
spiral structure is dominated by m = 2 spiral modes - we shall 
return to the nature of their angular momentum transport in the 
following sections. 

Figure[3]shows azimuthally averaged radial profiles for Simu- 
lations 1, 2, 3, 4 and 5. The impact of increasing angular momen- 
tum is evident in the surface density profiles (the top left panel of 
Figure[3ll. 

There appears to be a clear difference in trend between non- 
fragmenting discs (Simulations 1 and 2) and fragmenting discs 
(Simulations 3,4 and 5). In the non-fragmenting case, the surface 
density profile is visibly steeper (with a corresponding higher as- 
pect ratio, seen in the top right panel). As these plots are time- 
averaged, we can use each snapshot to derive the standard devi- 
ations of each curve for all radii, and perform fitting. All the 
curves are fitted to a pow er-law with exponential decay, of the form 
( IWilliams&Ciezall20Ilh : 



(where p and Rc are free parameters). This yields p ~ 1.99 for 
Simulation 1 and 2, and p ~ 1 for the other simulations. These re- 
sults are reflected in the Toomre Q profiles in the bottom left panel, 
where Simulation 2 displays a very narrow range of gravitation- 
ally unstable radii in comparison to the other discs. Indeed, despite 
their apparent similarities in the rest of Figure[3] Simulations 3, 4 
and 5 become distinct in their Q-profiles (although 3 and 4 remain 
similar). Simulation 5 displays the largest region of instability, ex- 
tending over a range spanning approximately 80 au. 

This set of results is an excellent example of disc fragmen- 
tation requiring at least two criteria be met. The bottom right 
panel shows the dimensionless cooling time parameter /3c, which 
again displays differing trends between fragmenting and non- 
fragmenting discs. Values of /3c of order unity present amenable 
conditions for fragmentation, and all four discs satisfy this condi- 
tion in their outer edges. Simulation 2 satisfies this at the lowest 
radii, at around 40 au. However, we must also note that at 40 au, 
Q ~ 3 and is increasing with radius. Therefore, while this region 
can maintain very cool gas, it is not gravitationally unstable, and 
we therefore do not expect it to fragment - this expectation is satis- 
fied, despite running the simulation for 30,000 years. As a test, we 
ran the disc for an extra 30,000 years, and no fragmentation was 
seen during this subsequent period. 

The observant reader will note that fragments are able to form 
in Simulations 3 and 4 at large radii, despite high values of Q. Fig- 
ure |4] goes some way to explaining this phenomenon. These discs 
are highly variable (Simulation 2 being the most steady-state of 
the four). While temperature fluctuations in these discs are rela- 
tively low, the fluctuations in Q are quite high, as has been observed 
previously in discs with infall ( iHarsono et alj|20Tl1) . Simulation 1 
maintains fluctuations of order 50% in the outer regions, ensuring 
Q » 2. The other simulations have extremely high Q fluctua- 
tions (nearly 100%), and all possess very low values of /3c, of order 
unity or less. It is this transient behaviour that allows Q < 2 for a 
sufficient period of time (say one orbital period at a given radius) 
which allows the disc to fragment. Simulation 1 also has extremely 
variable Q, but this translates into relatively weak surface density 
perturbations in regions of low E (as can be seen in the left panel 
of Figure|4]l, and therefore does not fragment. 

In effect, for Simulations 3 and 4 to fragment, the square 
of the epicyclic frequency must become negative, a phe- 
nomenon that typically defines the outer edge of accretion 
discs, where perturbations from circular flow will be unstable. 
This tends to agree with the consensus view that fr agmenta- 
tion i n the inner regions of discs is difficult to achieve jRafikovl 
'2005"; 'Matzner & Levin 2005; Whitworth & Stamatellos' l2006l: 
Bole v et al., ,2006: Stamatellos & Whitworth 2008 ; Forgan et al.1 
2009i; IClarkell2009l) . On the other hand, Simulation 5 shows less 
variation than the others, and satisfies the canonical fragmentation 
criteria (see next section) at around 100 au, where fragments are 
subsequently seen. The initial fragment provides an extra torque to 
the spiral wave that spawned it, which results in consecutive frag- 
mentation along the wave. For example, the condensation in the 
lower right quadrant of the lower right panel of Figure [2] subse- 
quently fragments into several objects. In all the simulations con- 
ducted, if multiple fragmentation occurs, it seems to occur as a re- 
sult of an initial fragment's torques on the disc, with time delays 
between the first object and the second object of up to an orbital 
period. 
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Figure 2. Surface density plots of the M^i = IMq simulations (Simulations 1, 2, 3, 4, 5 and 6). The snapshots were taken two free fall times after the 
collapse was initiated (for Simulations 1 and 2, top), or at the point of disc fragmentation (for Simulations 3, 4, 5 and 6, middle left and right, and bottom left 
and bottom right respectively). 



3.2.2 Angular Momentum Transport 

We have seen that the spiral structure appears to be dominated by 
m = 2 modes, and that the discs are highly variable. T hese are both 
symp toms of non-local angular momentum transport jporgan et al.l 
uOlli) . an d have been seen in previous simulation s with envelope 
accretion jKratter et al]|2008l : iHarsono et al]|201lh . We carried out 
a Fourier decomposition to study the amplitude of the m-modes. 
Figure |5] shows spectrally averaged m-modes as a function of disc 
position for Simulations 1-5. These confirm our visual inspection 



of the discs' spiral structure - the dominant m-mode in the unstable 
region of all five discs is m = 2, with the majority of the power 
remaining in the lower modes at larger radii. However, Simula- 
tion 1 shows moderate power across a wide range of mode num- 
berfl suggesting that it may not be subject to non-local transport. 



Spectral averaging weights the amplitude of each mode by m, giving rise 
to seemingly high average m modes as shown here - it is likely that m « 4 
modes dominate in the case of Simulation 1. 
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Figure 3. Azimuthally averaged radial profiles from the Afci = IA/q simulations (Simulation 1 (solid line). Simulation 2 (dotted lines). Simulation 3 (dashed 
lines). Simulation 4 (dot-dashed lines) and Simulation 5 (triple-dot-dashed lines)). Simulation 6 is excluded as it forms a multiple system. The profiles are 
time averaged from the formation of the disc to the point of fragmentation (in the case of Simulations 1 and 2, they are time averaged for approximately 7000 
years, coiresponding to approximately 19 ORPs). The top left panel shows the surface density profile (calculated from the material contained within the first 
three scale heights only), the top right shows the aspect ratio, the bottom left shows the Toomre Q parameter, and the right hand panel shows /3c. 




Figure 4. Variation in the mean surface density (left) and Toomre parameter Profile (right), for the Md = I.OMq simulations (1, 2, 3, 4 and 5). 



T his is confirmed by calculating the non -local transport fraction 
^ jCossins et al1l2009l : iForgan et alj|20l"ll) . All simulations except 
Simulation 1 exhibit values above 0.5, indicating non-local trans- 
port is significant. By contrast, Simulation 1 maintains values less 
than 0.2, indicating the transport is indeed local. 

Figure |6] shows azimuthally averaged profiles of the a- 
parameter, calculated for various components of the disc stress. The 
solid lines indicate the total disc stress, the dashed lines show the 
contribution from gravitational stresses Ofgrav, and the dotted lines 
show the Reyn olds stress QReyn (expressions for these functions 
can be found in lPorgan et al.l201ll) . We also plot the stress induced 
by the artificial viscosity, Oart to confirm that it is small in compar- 
ison to the other disc stresses. 

In contrast to isolated disc simulations, we can see that the 



Reynolds stresses are typically large, of order q ~ L This 

is in keeping with values obtained by iKratter et in ( l2008h and 
lHarsonoetallj201lh in similar disc configurations. 

Again, there is a stark contrast between non-fragmenting and 
fragmenting discs. We can see that the standard Qtot = agrav + 
ctRcyn > ctcrit =0.06 Criterion (derived for isolated discs) is in- 
valid. The large Reynolds stresses (induced by velocity shear from 
infalling envelope material contacting the disc's rotating flow) pre- 
vent the gravitational stress from dominating unless it can reach 
values closer to unity (in some cases up to 5 times larger than cvcrit)- 
We should therefore be careful about using this criterion in young 
systems that are still surrounded by a substantial envelope - indeed, 
we should develop ne w criteria that take this accretion into account 
Jporgan & Ric3l201lh . 
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Figure 5. Spectrally averaged m modes as a function of radius for the 
A4i = 1.0A/q simulations (1, 2, 3, 4 and 5). Note that inside 10-20 au 
the discs are dominated by numerical viscosity, meaning that data in this 
region should be ignored. With the exception of Simulation 1 , all discs ex- 
hibit non-local angular momentum transport. Consequently, the m = 2 
mode dominates consistently at large radii for Simulations 2-5. 



From this data, we suggest that fragmenting discs require the 
local gravitational stress to be of the order (or larger than) the local 
Reynolds stress. We can see that Simulation 2 satisfies this con- 
dition weakly at 40 au, again where Q is too large for the disc to 
be sufficiently unstable. At larger radii, the gravitational stress de- 
creases with distance, further ensuring no fragmentation can occur. 
Simulation 1 appears to satisfy this at 100 au, but ctgrav varies quite 
strongly, only dominating at much larger radii, where the surface 
density becomes negligible. The differences between Simulations 
1 and 2 are due to the different modes of angular momentum trans- 
port at play (local and non-local respectively). 

Simulations 3 and 4 show gravitational stresses comparable 
to local Reynolds stresses at around 150 and 110 au respectively 
- these correspond well to the initial formation radius of the frag- 
ments. Simulation 5 appears to be qualitatively different, in that it 
fragments in a regime where the gravitational stress is still small 
compared to the Reynolds stress, but Ograv > Qcrit- The disc ex- 
hibits more power in higher m-modes than the others, and the total 
Q-profile differs also. Where Simulations 3 and 4 show relatively 
flat atot outside 50 au. Simul ation 5 shows a stea dily increasing 
curve. Naively comparing with lForgan et al.l ( 1201 ih 's study of iso- 
lated discs, this would suggest that Simulation 5 shows more evi- 
dence of local angular momentum transport than Simulations 3 and 
4. However, as has been said. Simulations 3-5 all possess large non- 
local transport fractions. 

Simulations 3-5 produce disc fragments at a minimum dis- 
tance of 100 au from the central star, where the local (time- 
averaged) agrav at the fragmentation radius ranges from 0.06 to 
0.7. The subsequent evolution of the fragments in the non-linear 
regime vary between each simulation - for example. Simulations 
3 and 5 produce several subsequent fragments each, as a result of 
torques induced by the presence of the initial fragment, which itself 
undergoes rapid inward migration. Simulation 4 produces a single 
fragment, which grows to form a second pointmass, where the mass 
ratio M2/M1 ~ 0.15. 



4 DISCUSSION & CONCLUSION 

From a set of numerical experiments with simplified initial condi- 
tions, we have shown that the initial cloud angular momentum tends 
to be the critical factor in subsequent disc fragmentation (com- 
pared to t he initial thermal energy of the cloud). We confirm the 
results of iRice et al. I ( I2OIOI) , in that fragmentation does not occur 
if A < 10"''. While they also state that fragmentation occurs if 
A > 10"'', their parameter study does not resolve the fragmenta- 
tion boundary to the same level that this paper does. In this sense, 
we can say that their conclusions are broadly correct. Agreement 
can also be found on the fragmentation radius, which we also see 
increasing with initial cloud mass, and the steep surface density 
profiles obtained for non-fragmenting discs, placing a large amount 
of material within 20 au (although on this second statement we 
must caution that numerical effects dominate within 10 au). This 
agreement is somewhat surprising, in that they assume that the an- 
gular momentum transport due to self-gravity is local, where it has 
been established in this paper and i n others that it is lik ely to be 
dominated by global spiral modes (cf lHarsono et alfcoi ll) . 



All fragments form at distances above 70 au, typically con- 
sistent w ith pr evious studies of both is olated and accreting discs 



dRafikov 


I2OO5I; iMatzner & Levin '2005; Whitworth & StamatellosI 


20061; 


Bolevetal. 20061; Stamatellos & Whitworth ,2o6a 


Forean et al. 20091; Iciarke |2009|; Rice et alj 2010l). which in- 



dicate that, despite maintaining low values of Q, the inner regions 
do not cool efficiently or iinpart sufficient stress due to the 
gravitational instability to produce fragments. This is true in spite 
of non-local angular momentum transport or accretion from the 
envelope. It appears that maintaining surface density profiles 
which place large amounts of (cool) material at large enough 
radii is the key to fragmentation. The disc-to-star mass ratios are 
typ ically of order 0.5 at fra gmentation, although they can be lower 
(cf IStamatellos et alj|2oT ij). The masses of fragments formed are 
typically range from 1 to 15 Jupiter masses: these are roughly 
consistent with the local Jeans mass in the spiral waves from which 
they form (Forgan & Rice 2011), except in some cases where the 
Jeans length is sufficiently long to allow the break up of fragments 
into multiple objects. Future work will focus on the survival of 
such fragments as they condense out of the disc, to investigate 
the r ole of internal pressure support (^ Kratter & Murrav-Clavl 
I2OIII) and the importance of fragment-fragment collisions 
dShlosman & Begelmaj|l989h . 

Marginally unstable, self-gravitating discs which do not frag- 
ment can remain massive for long periods of time. In one case, the 
disc was shown to remain at masses roughly half that of the par- 
ent star for around 50,000 years after formation, showing that the 
self-gravitating phase could be more long-lasting than previously 
thought, much as the Class and Class I timescales have also been 
shown to be much longer than expected (,Evan&.201 li) . 

The inherent variability seen in all simulations conducted in- 
dicates that accretion onto the central star varies significantly. Even 
in the more quiescent, non-fragmenting simulations, the value of Q 
varies by up to 50% or more, ensuring episodes of gravitational in- 
stability marked by high-amplitude, low-m spiral waves followed 
by episodes of relative stability with weak spiral structures. This 
variability occurs despite the omission of stellar radiative feedback 
physics, unlike e.g. Stamatellos et al. ( 201 1), who demonstrate that 
accretion luminosity can enhance this episodic nature of instability. 
This would suggest that self-gra vity alone may account for the vari- 
abilit y seen in Class I protostars jBolev & Durisenl2008l ; l\^robvovl 
I2OO9I) . It is quite possible that the one-parameter boundary we find 
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Figure 6. Azimuthally averaged o-parameter for the Md = IMq simulations (Simulation 1 (top left panel). Simulation 2 (top right panel). Simulation 3 
(bottom left panel) and Simulation 4 (bottom right panel)). Simulation 5 is excluded as it forms a multiple system. The profiles are time averaged from the 
formation of the disc to the point of fragmentation (in the case of Simulation 1, it is time averaged for 7000 years, coiTesponding to approximately 19 ORPs). 
The solid line shows the total a, the dashed lines show agrav, the dotted lines show aji^y„, and the dot-dashed lines show the contribution to the stress from 
artificial viscosity, aart- Note that this numerical a dominates within around 10 AU, so we must disregard data from within this region. 



disappears with the addition of stellar luminosity, as its strength 
will clearly depend on Mc\. 

Indeed, this simple boundary is unlikely to apply in nature due 
to other complicating factors, including non-uniform density pro- 
files and the addition of turbulence. A density profile that decays 
with radius, such as a supercritical Bonnor-Ebert sphere is a more 
physical choice, ensuring more mass will be concentrated in the in- 
ner disc at early times. We suggest that these steeper profiles will 
require an increased amount of initial angular momentum to pro- 
duce fragmentation, and that the boundary established for uniform 
clouds is a lower limifl 

The addition of turbulence tends to weaken any correlation 
with angular momentum ( Walch et al. 2010). Spherical symmetry 
is broken by construction, and discs accrete from the environment 

* Strictly, an even lower limit could exist if the density increased with ra- 
dius, but these conditions are even less likely to appear in nature than simple 
uniform conditions 



through filamentary structures, allowing significant specific angu- 
lar momentum to be injected asymmetrically, exciting low-m spiral 
modes as the disc attempts to redistribute itself. In general, discs 
formed from turbulent simulations tend to grow more slowly, and 
have smaller radii for a given angular momentum - this again sug- 
gests that the criterion derived in this paper remains a lower limit. 



Considering that the threshold angular momentum required 
(A > 5 X 10~^) is sma ll compared to observa tions, which indicate 
a median of A ~ 0.02 jGoodman et"ai]|l993h - this would suggest 
that disc fragmentation will be a somewhat common process, but 
occur o n such short timescales t hat observations of such events will 
be rare jStamatellos et alj|201lh . Exactly how this threshold might 
vary with surface density profile and the addition of turbulence is 
difficult to quantify, but merits further research. 
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